capture program drop GenPolynomial
program define GenPolynomial
gen above           = totmuj>=20
gen totmuj2         = (totmuj-20)^2
gen totmuj3         = (totmuj-20)^3
gen totmuj_above    = (totmuj-20)*above
gen totmuj2_above   = (totmuj-20)^2*above
gen totmuj3_above   = (totmuj-20)^3*above
gen totmuj_original = totmuj
replace totmuj      = totmuj-20
end

graph drop _all 
use DTA/Datos_ENIA_95_07_Processed.dta, clear
GenPolynomial

global graph_vars km_total mw kw_total 

global Subsamples all period3 sectorB large 

local range = 10

foreach subsample in $Subsamples {
 preserve
 keep if abs(totmuj) <= `range' & `subsample'==1
   /*expand 2 if totmuj==0, gen(duplicateind)
> gen SE_obs = 1 if duplicateind==1
> replace totmuj          = -0.0001 if duplicateind==1
> replace totmuj_original = 20 - 0.0001 if duplicateind==1
> replace above           = 0 if duplicateind==1
> replace totmuj2         = (totmuj)^2 if duplicateind==1
> replace totmuj_above    = 0 if duplicateind==1
> replace totmuj2_above   = 0 if duplicateind==1
> replace totmuj_original = 20-0.0001 if duplicateind==1*/
 gen SE_obs=.
   foreach var of varlist $graph_vars {
   reg `var' totmuj totmuj2 totmuj_above totmuj2_above above i.ciiu3 i.year if SE_obs==., vce(robust)
   qui estat vce
   matrix B=r(V)
   matrix B=B[1..5,1..5]
  gen pl_`var'=totmuj*_b[totmuj]+totmuj2*_b[totmuj2]+totmuj_above*_b[totmuj_above]+totmuj2_above*_b[totmuj2_above]+above*_b[above]
  sum `var' if totmuj>=-3 & totmuj<0
  replace pl_`var'=r(mean)+pl_`var'
  gen pl_sd_`var'=totmuj^2*B[1,1]+totmuj2^2*B[2,2]+totmuj_above^2*B[3,3]+totmuj2_above^2*B[4,4]+above^2*B[5,5] ///
                                 +2*totmuj*totmuj2*B[1,2]+2*totmuj*totmuj_above*B[1,3]+2*totmuj*totmuj2_above*B[1,4]+2*totmuj*above*B[1,5] ///
                                 +2*totmuj2*totmuj_above*B[2,3]+2*totmuj2*totmuj2_above*B[2,4]+2*totmuj2*above*B[2,5] ///
                                 +2*totmuj_above*totmuj2_above*B[3,4]+2*totmuj_above*above*B[3,5] ///
                                 +2*totmuj2_above*above*B[4,5]
 replace pl_sd_`var'=sqrt(pl_sd_`var')
 }
  collapse $graph_vars pl_*, by(totmuj_original) // they have no variation anyway, only the actual data does
 foreach var of varlist $graph_vars {
 cap gen pl_sd_`var'_L = pl_`var'-1.96*pl_sd_`var'
 cap gen pl_sd_`var'_H = pl_`var'+1.96*pl_sd_`var'
 twoway  (scatter `var' totmuj_original , legend(order(1 "Real data" 2 "Predicted Polynomial" 4 "Prediction S.E.")) xtitle("Number of Women") ytitle("log-`var'") ) || ///
                 (line pl_`var' totmuj_original if totmuj_original<20, xline(20 ) lcolor(black) ) || ///
                (line pl_`var' totmuj_original if totmuj_original>=20, xline(20) lcolor(black) ) || ///
                 (line pl_sd_`var'_L pl_sd_`var'_H totmuj_original if totmuj_original<20, lpattern(dash dash) lcolor(gray gray)) || ///
                 (line pl_sd_`var'_L pl_sd_`var'_H totmuj_original if totmuj_original>=20, lpattern(dash dash) lcolor(gray gray)) ///
                 , name(G`var'_range`range'_subs`subsample', replace) graphregion(color(white))
 graph export GRAPHS/RDD_`var'_range`range'_sample`subsample'.pdf, replace
  }
  restore
  }